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The finite element ablation and thermal response (FEAtR, hence forth called FEAR) 
design and analysis program simulates the one, two, or three-dimensional ablation, internal 
heat conduction, thermal decomposition, and pyrolysis gas flow of thermal protection system 
materials. As part of a code validation study, two-dimensional axisymmetric results from 
FEAR are compared to thermal response data obtained from an arc-jet stagnation test in 
this paper. The results from FEAR are also compared to the two-dimensional axisymmetric 
computations from the two-dimensional implicit thermal response and ablation program 
under the same arcjet conditions. The ablating material being used in this arcjet test is 
phenolic impregnated carbon ablator with an LI-2200 insulator as backup material. The 
test is performed at the NASA, Ames Research Center Interaction Heating Facility. 
Spatially distributed computational fluid dynamics solutions for the flow field around the 
test article are used for the surface boundary conditions. 


Nomenclature 


2-D 

= 

two-dimensional 

3-D 

= 

three-dimensional 

Bt 

= 

pre-exponential factor for the i th resin component 

c p 

= 

solid material specific heat, J/kg-K 

Ch 

= 

Stanton number for heat transfer 

CFD 

= 

computational fluid dynamics 

DPLR 

= 

data parallel line relaxation 

F 

^ ai 

= 

activation energy for the i th resin component, J/kg-mole 

FEAR 

= 

finite element ablation and thermal response design and analysis program 

FIAT 

= 

fully implicit ablation and thermal response program 

h g 

= 

enthalpy of pyrolysis gas, J/kg 

h 

= 

weighted average solid enthalpy, J/kg 

H 

= 

thermocouple depth, m 

i 

= 

node index, resin component index (A, B ,C) 

IHF 

= 

interaction heating facility 

k 

= 

orthotropic thermal conductivity matrix, W/m-K 

K 

= 

pyrolysis gas mass flux vector, kg/m 2 -s 

MSL 

= 

Mars Science Laboratory 
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interpolation, or shape functions 
unit vector normal to the surface 
pressure, N/m 2 

phenolic impregnated carbon ablator 
net heat flux to the surface 
universal gas constant, J/kg-mole-°K 
thermocouple radial location, m 

recession rate vector, m/s 
temperature, °C 
time, sec 

two-dimensional implicit thermal response and ablation program 

thermal protection system 

boundary layer edge velocity, m/s 

weighting function 

parameter defined over one element 

absoptivity 

temperature difference, °C 
resin volume fraction 
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= permeability tensor 
= density of first resin component, kg/m 3 
= density of second resin component, kg/m 3 
= density of fiber reinforcement, kg/m 3 
= residual density, kg/m 3 

= solid material density, instantaneous density, kg/m 3 
= initial density, kg/m 3 
= density of pure virgin material, kg/m 3 
= density of pure char material, kg/m 3 
= density exponent factor 
= kinematic viscosity, m 2 /s 
= porosity 
= volume domain 
= gradient operator 
= denotes flux quantity, 1/m 2 


I. Introduction 

T HERE are two types of thermal protection systems (TPS) in use today, reusable and ablative. Use of reusable 
TPS, like the space shuttle tiles or a metallic heat sink, are generally limited to low heat flux atmospheric entry, 
hence forth called entry, trajectories. 1,2 Ablative TPS, on the other hand, can withstand large heat fluxes and heat 
loads, and are generally used for vehicles which have a high entry velocity, such as ballistic planetary entry 
missions. In addition to high entry velocity, ablative TPS are also well suited where the target planet has a high 
atmospheric density, such as Jupiter, since the heating rate is proportional to the density. 3 

Arcjet test facilities, such as the NASA, Ames Research Center Interaction Heating Facility (IHF), provide the 
primary means to test the performance of TPS materials 4 . Data obtained from these arcjet tests also provide a means 
with which to validate ablative thermal response computational tools and develop material response models. The 
finite element ablation and thermal response (FEAtR hence forth called FEAR) design and analysis program was 
developed to analyze ablative heatshields in multiple dimensions. 5,6 The one-dimensional validation and verification 
of FEAR was accomplished in Ref. 5. Verification of the three-dimensional (3-D) FEAR results was accomplished 
in Ref. 4 as well, but the lack of anisotropic material properties, computational fluid dynamics (CFD) solutions, and 
3-D arcjet data made full validation problematic, so only partial validation was achieved in Ref 5. Complete 
validation of the three-dimensional solution is the subject of future work. 


2 

American Institute of Aeronautics and Astronautics 



The ablating material used in this test was the phenolic impregnated carbon ablator (PICA), which is a low 
density TPS material developed in the 1990’s. 7,8 PICA was successfully flown on the Stardust sample-return 
capsule as the fore-body heatshield, and is the forebody heatshield for the Mars Science Laboratory entry vehicle. 
PICA was also one of the final two candidate TPS materials for the Orion Crew Exploration Vehicle heatshield. 
The Orion TPS Advanced Development Project conducted a large number of arc -jet tests over a wide range of 
conditions to evaluate the performance of PICA and validate its thermal response model. 

FEAR currently supports three and six noded triangle elements as well as four, eight, and nine noded 
quadrilateral elements. For this study, only four noded quadrilateral elements are compared to the arcjet test data 
and TITAN. The arcjet data includes measurements of surface temperature, surface recession, and in-depth 
temperatures. The heat flux, recovery enthalpy, and pressure boundary conditions are provided in the form a two- 
dimensional (2-D) spatial distribution from the data-parallel line relaxation (DPLR) CFD solution and is decoupled 
from the thermal response simulation albeit at the expense of losing some fidelity as the shape of the test article 
changes. In this paper, the 2-D FEAR solution is compared to both PICA arcjet data and solutions from the 2-D 
implicit thermal-response and ablation (TITAN) program 9,10 . TITAN is a 2-D, implicit finite difference solution 
scheme for the energy equation and three-component decomposition model with a moving grid. In this paper, the 2- 
D FEAR solution is compared to both PICA arcjet data and solutions from the TITAN program. 


II. Finite Element Ablation and Thermal Response Design and Analysis Program 

The governing differential equations which form the basis of the FEAR program describe the ablative thermal 
response problem and consist of the conservation of mass, momentum, and energy. A detailed derivation of the 
governing differential equations from first principles is presented in Ref. 5 and summarized in this section. The 
conservation of mass may be written as 

^- = -Vrh" g (1) 

dt g 

The instantaneous density is assumed to be a three-component mixture consisting of two resins and one reinforcing 
material components and is given by 


P= r (PA +P B )+( 1 ~ r )Pc 

The decomposition over time is calculated using an Arrhenius relation 11 ' 12 written as 


dt 



( 


pj-pj, 

Po, 


Vi 


i = A, B,C 


( 2 ) 


(3) 


In general, ablative materials are porous materials and the porosity, as well as the permeability, increase as the 
material decomposes from its virgin state to its charred state. Therefore, Darcy’s law is used to calculate the 
pressure drop due to the flow of pyrolysis gas through the material. 


m" vp 


V(f) 


(4) 


The conservation of energy for the conductive transport of energy through an ablative material that is 
decomposing in-depth may be written as 


pc n — =kVT+(h e -h)^+ m"Vh + S pcVT 

H p dt v g ’ dT g g p 


(5) 


III. Weak Formulation 

The corresponding weak form of the governing system of equations, (1), (4), and (5) may be constructed by first 
multiplying the residual of the differential equation by an appropriate set of weighting functions W and integrating 
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by parts over the domain Q. If the weighting functions W are set equal to the shape functions N which approximate 
the unknown variables, the result is known as the Bubnov-Galerkin form 13 . The method of weighted residuals 
requires that the weighted average of the error, or residual, vanish over the solution domain. The weighted residual 
statement for the conservation of mass and momentum is given in Eqs. (6) and (7), respectively. The weighted 
residual statement for the conservation of energy is given in Eq. (8). 


JIL 


iN^dQ. = 0 


dp 


dt 


E 


K„ 


m" (e) + -L vp (e) \N u ' v ' w dD. = 0 

8 V(j) 


... . . rfF (e) 

Jj| kVT {e) +(h g -h)-^-+m" g Vh g (e) +Spc p VT <e> -pc p — |y.JQ =0 




dt 


( 6 ) 

(7) 

( 8 ) 


In Eqs. (6) and (7), the shape functions that appear in the mass and momentum equations are different. In this 
formulation, the shape functions which approximate the pressure must be one order lower than those used to 
approximate the pyrolysis gas flux components in order to prevent an over constrained system of discrete equations. 
This formulation makes the pressure discontinuous across element boundaries and thus, pressure is not a primary 
variable in the weak formulation 14 . 

Integrating Eqs. (7) and (8) by parts to eliminate the pressure derivative and reduce the derivative on the 
temperature to first order results in Galerkin’s weak form for the conservation of mass, momentum and energy given 
by Eqs. (9), (10), and (11), respectively. 

Q. 

K* r)N u,v,w 
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(ii) 




p(e) 


The surface integrals introduced by the integration by parts provide the connection to the surface boundary 
conditions present on any surface of the element. In the momentum equation, the first surface integral represents 
surface tractions and the second a specified mass flux at the surface. In the energy equation, the surface integral 
represents the net conductive flux into the volume, from any surface bounded by T and can be expanded to 
represent any number of boundary conditions by summation over all surfaces bounded by T. Since Eq. (11) is for 
one element, T includes all of the surfaces of the element. For example, in the case of an ablative material, the 
surface boundary conditions are given by 

|| (q-n)N i dr = p e U e C H (H si -h sw + B’h c +B’ g h g -B\)-q* + q rad -aq rad ( 12 ) 

w» out in 


The governing differential equations for linear elasticity are used to solve for the thermal stress and also to move 
the mesh while recession is occurring. These governing differential equations and corresponding weak form have 
been developed and described in detail in several references and are not presented in this paper 13,15 . 
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IV. Arcjet Coupon Geometry and Computational Mesh 

One advantage in the FEAR program is that the geometry of the object being analyzed can be created in 
ProEngineer®, or any other computer-aided design (CAD) software which has the ability to export the geometry in 
IGES or STEP file format. Once the geometry is created, the geometry can be imported into MSC PATRAN for 
meshing, material specification, and boundary condition application. This portability facilitates the design process 
by allowing designers, and both structural and thermal analysts to have the ability to work with a common model. 
Data exchange is simplified since the FEAR solution can either be imported back into PATRAN or into Tecplot. 
The computational geometry being used in this paper was developed in ProEngineer and is based on arcjet test 
coupons designed and fabricated by NASA, Ames Research Center. 16 

The arcjet coupon is an Iso-Q shape with a 10.16 cm diameter and overall length of 9.78 cm. The coupon is 
fabricated from two materials, the test material which is PICA, and the backup insulation layer which is LI-2200, a 
low density silica fiber based ceramic insulator. 17 The PICA cylinder is 4.127 cm diameter in the middle section and 
5.08 cm diameter near the edges. The axisymmetric geometry, after being imported into PATRAN, is shown in Fig. 
1 with the computational mesh superimposed. The actual arcjet specimen has a 1.78 mm air gap between the PICA 
and the LI-2200 in the 4.127 cm depth section of the PICA, but the air gap is neglected in this analysis since surface 
to surface radiation is not currently supported in PEAR. After the geometry is imported into PATRAN, the 3-D 
solid is broken down into a cross-sectional set of surfaces and only those on one side of the centerline are retained to 
perform the axisymmetric analysis. 

The arcjet coupon geometry was meshed with the meshing routine built within MSC PATRAN®. The surfaces 
were meshed with structured, 4-noded quadrilateral elements. The mesh was biased towards the heated surface and 
consisted of 1188 nodes and 1119 elements. In the region near the shoulder of the coupon, the mesher produced 
much smaller elements in order to maintain congruency between the individual geometric surfaces which make up 
the axisymmetric geometry. Puture work will examine different element types such as triangles, wedges, and 
unstructured quadrilaterals, and also examine the performance of higher order elements. 


9.78 cm 



z X 


Fig. 1 PICA arcjet coupon geometry and computational mesh. 

The arcjet coupon had six thermocouples embedded in-depth in a 2.54 cm plug which was inserted down the 
centerline. The coupon also had thermocouples at three radial locations at two depths equal to the two deepest 
thermocouples along the centerline. The thermocouple locations are summarized in Table 1. 


Table 1 Thermocouple locations 


Designation 

Depth, H 

Radial Distance, r 


(cm) 

(cm) 

TCI 

0.381 

0.0 

TC2 

0.762 

0.0 
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TC3 

TC4 

TC5 

TC6 

TC7 

TC8 

TC9 

TC10 


1.143 

1.524 

2.286 

3.048 

2.286 

2.286 

2.286 

3.048 


0.0 

0.0 

0.0 

0.0 

2.54 

3.81 

4.45 

4.45 


V. Boundary Conditions 

The surface boundary conditions were obtained from a CFD solution of the coupon in the arcjet flow. The 
centerline enthalpy for this test was 19.3 MJ/kg and the centerline heat flux was 255 W/cm 2 . The coupon was 
exposed to the flow for 40 seconds and allowed to cool. The cool down time for analysis purposes was 600 seconds. 
The heat flux varied across the surface of the coupon and also wrapped around the shoulder where there was non- 
trivial sidewall heating calculated. The distributed surface heat flux is shown in Fig. 2 and is applied to the external 
edges of the PICA and LI-2200. 

In addition to the heat flux boundary condition, FEAR requires two supplementary boundary conditions. The 
first is a no flow boundary condition which is specified on the back edge of the PICA which interfaces with the LI- 
2200. The no flow boundary condition is required for the solution of the mass and momentum equations. The final 
required boundary conditions are to specify the zero displacement boundaries. The zero displacement boundaries 
are required for both the mesh movement scheme and the thermal stress calculation and need not be coincident. The 
zero displacement boundary condition for the mesh movement scheme encompasses all of the nodes and elements 
representing the LI-2200 since only the PICA is allowed to recede. The zero displacement boundary for the thermal 
stress calculation is the bottom edge of the LI-2200. Radiation in or out of any edge, (or surface in 3-D) is not a 
hard requirement, but is applied to all the external edges of the coupon except for the bottom edge of the LI -2200. 
The initial temperature is 20°C and the radiation sink temperature is 21.1 °C. 



Fig. 2 CFD predicted surface heat flux distribution. 

VI. Results and Discussion 

FEAR was run in 2-D axisymmetric mode for the arcjet conditions and geometry described in Sections IV and 
V. The temperature distribution at 40 seconds when the arcjet flow is cutoff is shown in Fig. 3. The temperature on 
the side of the coupon has increased significantly compared to the temperature along the centerline at the same 
depth, a clear indication that the sidewall heating has a significant affect. 


6 

American Institute of Aeronautics and Astronautics 


0.08 - 


Temperature °C 


0.06 - 


>■ 


0.04 - 


0.02 - 


oj- 




--F- — 

— 

— 

— 



mimn l l 








= 

:=t= 1 











— 

— 


— 

L_ 

■ m llll — 








II r in« 






V ,1 II’ 








W ////////ffiiiir 








W 








W 7////MHI II' 








I ///. flffll 








1 mu 
















r mi, nun 








hi 








in 








lllllllUll 








/iniiiiiii 








niiiiififi 
















mini 








iminimri 









0.02 0.04 0.06 0.08 


■ 

■ 


X 


2400 

2300 

2200 

2100 

2000 

1900 

1800 

1700 

1600 

1500 

1400 

1300 

1200 

1100 

1000 

900 

800 

700 

600 

500 

400 

300 


Fig. 3 Temperature distribution at arcjet flow cutoff (40 sec). 


The temperature distribution after cooling down for 560 seconds is shown in Fig. 4 and is superimposed on the 
original outer mold line of the coupon to show the shape change. The coupon maintains the Iso-Q shape, as 
expected, and the side wall heating is not significant enough to cause recession in the radial direction. Another 
indication that the side wall heating effects are significant is to examine the density distribution. The final coupon 
char fraction distribution is shown in Fig. 5. The char profile indicates that the PICA has decomposed significantly 
near the side of the specimen compared to the same depth along the centerline. The boundary elements between the 
PICA and LI-2200 share common nodes and as a consequence, the contour plot looks as if there is a severe density 
gradient in those elements. Computationally, each node has the correct density applied when the element is 
assembled into the global matrix. Until contact boundary conditions are implemented within FEAR and the 
boundary elements have independent nodes, this output anomaly for density will have to be tolerated. 



Fig. 4 Temperature distribution after a 560 second cool down. 
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Fig. 5 Final char distribution. 


Comparison of the centerline thermocouple arcjet data with the results from both FEAR and TITAN are plotted 
Fig. 6. There was no pyrometer data available for comparison, but from Fig. 6 both TITAN and FEAR calculate 
similar surface temperatures. Recession measured on the centerline at the stagnation point was 4.06 mm and the 
stagnation point recession calculated by FEAR was 3.28 mm, a difference of 0.78 mm. The in-depth centerline 
temperatures calculated by FEAR compare well with both the TITAN solution and the arcjet data. Comparison of 
the temperatures along the radial lines at depths of 2.286 cm and 3.048 cm are plotted in Fig. 7 and Fig. 8 
respectively. The results from FEAR compares well to both the arcjet data and to the results from TITAN. The 
radial comparison indicates that the anisotropic material properties and the application of sidewall heating within 
FEAR are being implemented correctly. 
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Fig. 6 Centerline thermocouple temperature comparison between FEAR and TITAN 


The results from FEAR do not match those of TITAN perfectly. The inputs used by both FEAR and TITAN are 
referred to as a response model. The differences between FEAR and TITAN arise because of how the response 
model is developed. The thermal response model used in this study was taken from the Crew Exploration Vehicle 
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(CEV) Thermal Protection System (TPS) Advanced Development Project (ADP) material property database. This 
thermal response model was developed with measured material property data and numerous instrumented arcjet 
tests. The fully implicit ablation and thermal -response (FIAT) program, and the 2-D counterpart TITAN, benefit 
from the fact that adjustments to the thermal response model are made such that the FIAT solution matches the 
arcjet data over a wide range of conditions. Typically, the char thermal conductivity is adjusted during this process 
to account for all the other uncertainties present in the model, but other adjustments, such as altering the pyrolysis 
gas and virgin material compositions are also made. 



Fig. 7 Temperature comparison between FEAR and TITAN along a radial line at a depth of 2.286 cm. 



Fig. 8 Temperature comparison between FEAR and TITAN along a radial line at a depth of 3.048 cm. 

The physics being modeled within FEAR are nearly identical to that of FIAT and TITAN except for how the 
pyrolysis gas flow is treated. FEAR solves Darcy’s Law to determine both the magnitude and direction of the 
pyrolysis gas flow and the internal pressure. As a consequence, FEAR depends on the porosity of the material, 
permeability and pyrolysis gas viscosity. The database for these three properties is severely lacking and some 
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assumptions which are described in Ref. 5, must be made to arrive at a solution. Essentially, FEAR is introducing 
additional uncertainties and the thermal response model needs to be adjusted to account for them. If given the 
opportunity to make adjustments to the thermal response model, FEAR could be made to match the data as well as 
TITAN. With this new analysis capability, it is imperative that porosity, through-the-thickness and in-plane 
permeability, and pyrolysis gas viscosity be measured. Also, a proper database must be developed to reduce the 
uncertainty in these properties and a response model specific to FEAR must be developed. 

VII. Conclusion 

A 2-D axisymmetric analysis of an arcjet coupon constructed of PICA with an LI-2200 insulator was performed 
with the thermal response code FEAR. The solution results from FEAR compared well to both arcjet data and to the 
2-D TITAN solution results. The thermal response code FEAR also calculated the stagnation point recession to 
within 0.78 mm of the measured recession. The anisotropic property representation implemented within FEAR was 
shown to be correct. The thermal response model required by FEAR may need to be adjusted slightly from the 
response model developed for TITAN in order for the calculations to correlate more closely with arcjet test data. 
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